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Abstract 

In a recent paper an infinitely long memory model (the t- 

model) for the Euler equations was presented and analyzed. The model 
can be derived by keeping the zeroth order term in a Taylor expansion 
of the memory integrand in the Mori-Zwanzig formalism. We present 
here a collection of models for the Euler equations which are based 
also on the Mori-Zwanzig formalism. The models arise from a Taylor 
expansion of a different operator, the orthogonal dynamics evolution 
operator, which appears in the memory integrand. The zero, first and 
second order models are constructed and simulated numerically. The 
form of the nonlinearity in the Euler equations, the special properties 
of the projection operator used and the general properties of any pro- 
jection operator can be exploited to facilitate the recursive calculation 
of even higher order models. We use our models to compute the rate 
of energy decay for the Taylor-Green vortex problem. The results are 
in good agreement with the theoretical estimates. The energy decay 
appears to be organized in "waves" of activity, i.e. alternating periods 
of fast and slow decay. Our results corroborate the assumption in JHj , 
that the modeling of the 3D Euler equations by a few low wavenumber 
modes should include a long memory. 

1 Introduction 

The advent of very powerful computers enhances our ability to probe com- 
plex systems and reveal their dynamics. Yet, there are many problems 
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where the present computational capacity is not enough to fully resolve all 
their dynamic features. This situation poses a great challenge for numerical 
analysts. How can one extract meaningful information about a system's 
evolution when the direct simulation of such a system is out of reach? One 
approach is to use the qualitative information gathered by analytical, nu- 
merical and experimental studies as a guide for the formulation of reduced 
models which will hopefully reproduce the qualitative behavior observed so 
far, and moreover, reveal something new about the behavior of the system. 
Fluid turbulence is the archetypal example of a complex system and one of 
the main subjects of research for scientific computing. Even though a lot of 
information has been collected regarding the qualitative behavior of a real 
fluid (and of the associated equations of fluid flow) |21 13) E| 171 \T2\ ITS) 12U] . 
this information is very difficult to incorporate in the construction of re- 
duced models. In the present paper we construct a collection of reduced 
models for the 3D Euler equations which describe the evolution of an invis- 
cid fluid, based on the physical and numerical observation that the evolution 
of very smooth initial conditions can give rise to organized structures (known 
as vortices, vortex filaments, pancakes, sheets etc.) |15| . These organized 
structures are known to create long temporal correlations and since tem- 
poral correlations appear as memory integrands in our reduction formalism, 
the models we will construct are predominantly long memory ones. 

The reduced models we will construct are based on the Mori-Zwanzig 
formalism of irreversible statistical mechanics |221 1321 125| as reformulated 
by Chorin, Hald and Kupferman In this formalism, the correlations 
of the unresolved modes appear in the integrand of a memory term. In 
previous work it was found that the form of the memory term can 
simplify significantly if the approximation of a very long memory is made. 
This approximation can be effected by expanding in a Taylor series the 
memory integrand and keeping only the zeroth order term. That simplified 
model is known as the t-model and it was applied to the Euler equations jlUj 
and to the Burgers equation pj. Here we propose an alternative construction 
which is also based on a Taylor expansion of an operator (the orthogonal 
dynamics operator) which appears in the memory term. The t-model can 
be derived as special case of the present construction. We construct reduced 
models that retain the zeroth, first and second order of the Taylor expansion 
of the orthogonal dynamics operator. In addition, we derive a set of rules 
that can facilitate the recursive calculation of high order models. The rules 
are based on the observation that the form of the terms appearing in the 
reduced models is determined by the form of the nonlinearity, the special 
properties of the projection operator used and the general properties of any 
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projection operator. 

It is known that the formation of organized structures in a fluid flow is 
manifested through a cascade of energy from the large to the small scales 
(see e.g. ^5]). Even though the Euler equations conserve the energy of a 
smooth solution, they do not have to conserve the energy of a non-smooth 
solution. There exist estimates (see |3Uj and references therein) about the 
rate of decay of energy in an inviscid flow. We use our models to compute 
the rate of energy decay for the Taylor-Green vortex |311ll6j . The results are 
in good agreement with the theoretical estimates. It is interesting that the 
energy decay appears to be organized in "waves" of activity, i.e. alternating 
periods of fast and slow decay. 

The paper is organized as follows. In Section|2]we review the Euler equa- 
tions and the problem of underresolved computations for these equations. 
In Section |3] we give a brief presentation of the Mori-Zwanzig formalism 
which is the starting point of our approximations. In Section 0] we present 
the models based on the Taylor expansion of the memory term integrand. 
We also give a set of rules for the recursive evaluation of high order models. 
Numerical results of the models for the Taylor-Green vortex are presented 
in Section [5J Finally, in Section we discuss the numerical results and what 
they suggest for future work. 

2 The Euler equations and the problem of under- 
resolved computations 

Consider the 3D incompressible Euler equations with periodic boundary 
conditions in the cube [0,2-7r] 3 : 



vt + v-Vv = -Vp, V • v = 0, (1) 

where v(x,t) = (vi(xi,X2,X3,t),v 2 (xi,X2,x 3 ,t),v 3 (xi,X2,x 3 ,t)) is the ve- 
locity, p is the pressure and V = (gf^-, ^j)- The system in (^Q) is sup- 
plemented with the initial condition v(x,0) = vq(x) which is also periodic 
and incompressible and x = (x±, x%, x%). 

Since we are working with periodic boundary conditions, we expand the 
solution in Fourier series keeping M modes in each spatial direction, 

v M (x,t) = u k(t)e ikx , 

keFuG 
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where FU G = [- f , f - 1] x [- f , f - 1] x [- f , f - 1] . Also fc = (A* , A; 2 , k 3 ) 
and itfc(t) = (uj,(t) , uj,(t) , u^(t)) . We have written the set of Fourier modes as 
the union of two sets in anticipation of the construction of the reduced model 
comprising only of the modes in F = [-f,f-l]x [-y, ^ — 1] x [— ^, y — 1], 
where iV < M. The equation of motion for the Fourier mode becomes 

= _i fc ' u P A k u qi ( 2 ) 

where = / — is the incompressibility projection matrix and I is the 
3x3 identity matrix. The symbol • denotes inner product in M 3 . The system 
(J2J) is supplemented by the initial condition u$ = {ufc(O)} = {uofc}, k € 
FUG, where tto/c are the Fourier coefficients of the initial condition vq (x) . 

Even if we start from a very smooth initial condition and M is of the 
order of 10 3 in each direction (the state of the art in massively parallel 
computers), the solution of the system of ordinary differential equations (|2|) 
can create significant activity in the highest modes of our allowed resolu- 
tion. This phenomenon, called nonlinear instability, stems from the cascade 
(transfer) of energy from the large to the small scales and renders our cal- 
culations meaningless once a significant amount of energy has accumulated 
in the smallest scales of the solution. The source of the problem is that 
the system of equations @ conserves the energy E = \ ^2 |u/J 2 . Thus, 

keFUG 

when energy reaches the highest modes in our calculation, there is no way 
for this energy to exit the range of allowed wavenumbers, so it begins to 
contaminate the solution and leads to a catastrophic increase of the error. 
Hence, if one hopes to extract meaningful information from a finite calcula- 
tion, one needs to account for the required drain of energy out of the allowed 
range of wavenumbers (the same situation appears also in the case of the 
Navier-Stokes equations with small viscosity). Since we assume that we can 
only afford to solve the system (j2J) for M modes in each direction (the set 
FUG), our task is to construct a reduced model for the modes in F and use 
the modes in G to effect the needed drain of energy out of the set F. This 
is why we wrote the set of Fourier modes as a union of two sets. From now 
on, the set F will be called the resolved modes and the set G the unresolved 
modes. 

A lot of work (e.g. [E2 EJ ESI Ell EHl EHl) has been devoted to the 
development of models that effect the drain of energy needed to keep the 
calculation well resolved. The problem with such models is that they are 
introduced ad hoc, involve parameters that need to be adjusted and usually 



4 



work for some cases but not others. Another usual drawback is the per- 
turbative treatment of the nonlinear term. While this is adequate for large 
values of the viscosity, it is not for the more realistic small viscosity cases 
(see |29j and references therein). There is an obvious need for models that 
are derived directly from the Euler equations based on assumptions that re- 
spect the observed physics. In order to do that one needs as a starting point 
a framework that allows, in principle, the construction of reduced models of 
arbitrary precision. After such a formalism has been established, reduced 
models of varying precision can be constructed by incorporating qualitative 
knowledge about the behavior of the solutions. One instance of such a for- 
malism is the Mori-Zwanzig formalism of irreversible statistical mechanics 
[SH US] as reformulated by Chorin, Hald and Kupferman [SJIE1- As with 
every other formalism that allows the construction, in principle, of models 
of arbitrary accuracy, the actual formulation of the reduced model requires 
a lot of information about the full system. This could appear to defeat the 
purpose of constructing a reduced model. However, it is to be expected, be- 
cause each system has a certain amount of information and a reduced model 
should account for this information. In ^Hj, starting from the Mori-Zwanzig 
formalism and observations about the behavior of the solutions of fluids for 
very small or zero viscosity, an infinitely long memory model was proposed 
for the Euler equations. Such a model was based on the observation (see 
e.g. pQ) that the formation of vortices in a fluid flow can lead to very slowly 
decaying velocity temporal correlations. Since correlations are the building 
blocks of the memory term in the Mori-Zwanzig formalism, the assumption 
of an infinite long memory, which simplifies considerably the construction 
and the form of the reduced model, was natural. We should note here that 
the reduced model in (known as the t-model) goes against the common 
practice in modeling, where usually the opposite assumption of extremely 
short, or none at all, memory is invoked. As was shown in [11] . the two 
cases, of extremely long and extremely short memory, are the two sides of 
the same coin. However, there is a huge difference between them as far as 
their domain of validity. In the very short memory model, the unresolved 
modes are assumed to have very fast decaying temporal correlations, while 
the opposite is assumed in the case of the very long memory model. Such 
locality or non-locality of the correlations should ultimately be determined 
by the choice of resolved modes |27| . 
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3 The Mori-Zwanzig formalism 



We give a brief account of the Mori-Zwanzig formalism ([5]) which is the 
starting point for the different models to be presented in the next section. 
Note that here we use the standard notation associated with the Mori- 
Zwanzig formalism, and the use of the variables u, x should be clear from 
the context. 

Suppose we are given an M-dimensional system of ordinary differential 
equations 

with initial condition (f)(0) = x. 

The system of ordinary differential equations we are asked to solve can 
be transformed into the linear partial differential equation 

u t = Lu, u(x,0) = g(x) (4) 

where L = YaLi Ri( x )'ut : ana - * ne solution of ® is given by u(x,t) = 
g(4>(x,t)). Consider the following initial condition for the PDE 

g(x) = Xfc u(x,t) = 4> k {x,t) 

Using semigroup notation we can rewrite (jlj as 

d 

— e tL x k = Le tL x k 
at 

Suppose that the vector of initial conditions can be divided as x — 
where x is the iV-dimensional vector of the resolved variables and x is the 
(M — iV)-dimensional vector of the unresolved variables. Let P be an or- 
thogonal projection on the space of functions of x and Q = I — P. In previous 
publications on the Mori-Zwanzig formalism, the projection operator P was 
defined through a probability density function f(x) on the set of variables 
x. For the special case of the Euler equations, such a density is impossible to 
come by analytically or experimentally, thus we choose a different projection 
operator, which does not require the knowledge of a density for the values 
of x. For a function h(x) of all the variables, the projection operator we 
will use is defined by P(h(x,x)) = h(x,0), i.e. it replaces the value of the 
unresolved variables x in any function h(x) by zero. Similarly, the initial 
condition x = (x,x) is replaced by (x, 0). We should note here that such a 
projection operator is rather natural for the problem of the evolution of a 
very smooth initial condition, where one expects only a few Fourier modes 
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to have nonzero values initially. If we divide the wavenumbers in shells of 
different radii, then we can order the Fourier modes depending on which 
shell they belong to. For the numerical examples in this paper, the resolved 
modes will be taken in the first few shells, while the unresolved modes in 
the rest of the shells allowed by our resolution. 
The equation Q can be rewritten as 

^-e tL x k = e tL PLx k + e tQL QLx k + I e {t ~ s)L PLe sQL QLx k ds, (5) 
(Jt Jo 

where we have used Dyson's formula 

e tL = e tQL + /"* e (t-s)Lp Le sQL dSi (g) 

Jo 

Equation (jSJ is the Mori-Zwanzig identity. Note that this relation is exact 
and is an alternative way of writing the original PDE. It is the starting point 
of our approximations. Of course, we have one such equation for each of the 
resolved variables 4> k ,k = 1, . . . ,N. The first term in © is usually called 
Markovian since it depends only on the values of the variables at the current 
instant, the second is called "noise" and the third "memory". The meaning 
of the different terms appearing in (J5J) and a connection (and generalization) 
to the fluctuation-dissipation theorems of irreversible statistical mechanics 
can be found in [TT] . 
If we write 

e tQL QLx k = w k , 

w k (x,t) satisfies the equation 

■§iw k (x, t) = QLw k (x, t) 
w k (x,Q) = QLx k = R k {x) - PR k (x). 

If we project Q we get 

d 

P—w k (x,t) = PQLw k (x,t) = 0, 

since PQ = 0. Also for the initial condition 

Pw k (x,0) = PQLx k = 

by the same argument. Thus, the solution of (J7J) is at all times orthogonal 
to the range of P. We call (JJJ) the orthogonal dynamics equation. Since 
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the solutions of the orthogonal dynamics equation remain orthogonal to the 
range of P, we can project the Mori-Zwanzig equation (jSJ) and find 

^-Pe tL x k = Pe tL PLx k + P f e (t " s)L PLe sQL QLx k ds. (8) 
ot Jo 

4 Mori-Zwanzig models based on the Taylor ex- 
pansion of the orthogonal dynamics operator 

The t-model E3 can be derived in different ways but the one that is closer 
to our approach consists of expanding the memory integrand e^~ s ^ L PLe s ® L 
around s = and retaining only the zero order term. The essence of the 
t-model approximation is the approximation of the orthogonal dynamics op- 
erator e t( * L by the full dynamics operator e tL . For the models presented here, 
we proceed in an alternative way by expanding the orthogonal dynamics op- 
erator around s = 0. Depending on how many terms we keep (1,2,3, . . .), 
we obtain zeroth, first, second, .... order approximations respectively, 

s 2 

e s Q L = I + sQL + —QLQL + 0{s 3 ), (9) 
s 2 

PLe sQL = PL + sPLQL + —PLQLQL + 0(s 3 ). (10) 

Every term in the expansion has one more factor of QL than the previous 
term. At first sight, it seems that the construction of high order terms is 
just a tedious series of differentiations testing one's stamina. However, there 
are a few simple rules that can be used to facilitate the derivation of high 
order terms (see Section [4.2j) . 

There is no reason to expect a priori that the orthogonal dynamics oper- 
ator can be written as a Taylor series of s. In fact, we expect the dependence 
of e s ® L on s to be rather rough due to the rapid energy cascade to the unre- 
solved modes. The expression of this cascade is the fact that the orthogonal 
dynamics operator evolves quantities in a way that they acquire a compo- 
nent in the orthogonal complement of the range of the operator P. However, 
the results obtained from models of different orders in s are instructive and 
worth presenting (see also the discussion in Section EJ). 

Before we present the different models, we rewrite the equations © to 
conform with the Mori-Zwanzig formalism. We set 

Rk(u) = -i k ■ u p A k u q 

p+q=k 
p,q£FuG 
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and we have 

d -^ = R k (u) (li) 

for k G F U G. The system (|lip is supplemented by the initial condition 
uo = (uq,uq) = (no,0). Note that we focus on initial conditions where the 
unresolved Fourier modes are set to zero. This is enough, as mentioned 
before, for our purposes, since the evolution of even very smooth initial 
conditions can still give rise to the phenomenon of nonlinear instability. We 
will proceed and construct reduced equations for the Fourier modes u k with 
k G F. Of course, these equations will depend on the values of the Fourier 
modes in G and so, our task in constructing a model for the modes in F, 
is to model the behavior of the modes in G so that we can obtain a closed 
system for the modes in F. 



4.1 Zeroth, first and second order models 

We start our presentation of the models with the zeroth order model which 
is cubic in the Fourier modes. 

^Pe tL u 0k = Pe tL R k (u ) + P J* e^ L Z k (u )ds, (12) 

where 

Zl(u ) = PLQLu ok = -i( ^ k ■ Rp(u )A k u 0q + ^ k ■ u 0p A k R q (u ) J 

p+q=k p+q=k 
peG, q€F p€F, q€G 

(13) 

and 

Rk(u ) = R k (u ,0) = -i ^2 k ■ uo p A k uo q . 

p+q=k 
p,q£F 

If we make the approximation J** e^~ s ^ L Z^(uo)ds ~ f Q e tL Z%,(uo)ds = te tL Z k (uo), 
then we recover the t-model of jlUj . 

The system in (j!2[) is not closed in the quantities Pe tL UQ k for k G F. 
The reason is that the projection operator P does not in general commute 
with the nonlinear functions e tL R(u,0) and e tL Z°(u, 0). The simplest way 
to obtain a closed system is to commute the nonlinear functions and the 
projection P. This is a kind of mean-field approximation. The difference 
with other mean-field approximations is that here we do account for the 
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fluctuations of the unresolved modes, through the inclusion of the memory 
term. Thus, the reduced system becomes 

§- f Pe tL n k = Rk(Pe tL n ) + J* Z° k (Pe^ L u )ds, (14) 

The first order model is quartic in the Fourier modes and is given by 

^Pe tL u ok = R k (Pe tL u ) + jf Z° k (Pe^ L u )ds (15) 

+ f sZl(Pe^ L u )ds, (16) 
Jo 



where 



^fcC"o) = PLQLQLuo k = 
^ k ■ R p (uo)A k R q (u ) + ^2 k- Rp(u )A k R q (u )+ 

p+q=k p+q=k 

peFuG, qeG peG, qeFuG 

^2 k ■ Z°(u )A k u 0q + ^ k ■ uopA k Z°(u )\ 

p+q=k p+q=k 
peG, qeF p&F, q&G 



The second order model is quinitic in the Fourier modes and is given by 

^Pe tL u ok = R k (Pe tL u ) + £ Z° k (Pe^ L u )ds (17) 

+ f sZ 1 k (Pe^ L u )ds+ f ^Z 2 k {Pe (t -^ L u )ds, (18) 
Jo Jo 2 
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where 

Z%(uo) = PLQLQLQLu ok 



if Yl k ■ Z° p (u )A k R q (u ) + k ■ R p (u )A k Z°(u )+ (19) 

V ^ i — u ^ i — u 



p+q=k p+q=k 
peFUG, qeG peG, qeFUG 

Y k ■ B p (u )A k R q (u ) + Y k ■ R p (u )A k B q (u )+ (20) 

p+q=k p+q=k 

peFuG, qeG peG, qeFuG 

Y k ■ Z°(u )A k R q (u ) + Y k-R p (u )A k Z° q (u )+ (21) 

p+q=k p+q=k 

peG, qeF peF, qeG 

Y k ■ Z°(u )A k R q (u ) + Y k ■ R p (u )A k Z°(u )+ (22) 

p+q=k p+q=k 

peFuG, qeG peG, qeFuG 

Y k-Z°(u )A k R q (u )+ Y k-R p (u )A k Z°(u )+ (23) 

p+q=k p+q=k 

peG, qeFuG peFuG, qeG 



Y k ■ Z p (u )A k u q + Y k ■ u pA k Z q (u ) J 

k „ _i_ „—u / 



(24) 



p+q=k p+q=k 

peG, qeF peF, qeG 

where 

Bk(u ) = -i Y k ' R p{uo)A k UQ q . 

p+q=k 
p,qeF 

Some terms in p9 |) -(|24 |l can be grouped together to yield a simpler 
expression. Yet, it is more instructive to keep them separate because they 
reveal the rules that determine the types of expressions that appear in a 
term of arbitrary order. 

One common aspect of all the models of order one and higher, is that they 
involve integrodifferential equations where the integrals are convolutions. 
The computation of convolution integrals can be very expensive. However, 
the form of the convolution integral in our case can be translated into a sum 
of ordinary integrals by a simple change of variables. For the first order 
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model (fT3)) the change of variables s' = t — s gives 



^Pe tL u ok = R k (Pe tL u ) + [ Z° k (Pe sL u )ds + [ (t - s)Zl(Pe sL u )ds 

Jo Jo 



dt 



R k (Pe tL u ) + [ Z^Pe sL u )ds + t [ Zl(Pe sL u )ds - [ sZl(Pe sL u )ds, 
Jo Jo Jo 

The resulting integrals are no longer of convolution type and need not be 
evaluated from scratch for each value of t. Instead, they can be evaluated 
by adding the contribution of each new timestep to the existing value of 
the integral. Similar constructions can be used for the convolution integrals 
appearing in the higher order terms. 

Instead of expanding the orthogonal dynamics operator, one can start 
with equation (JSJ 

^-Pe tL u ok = Pe tL PLu ok + P [ e sL PLe {t ~ s)QL QLu ok ds. 
dt Jo 

(note the change of variables s' = t — s) and build an infinite hierarchy 
of Markovian equations which is equivalent to the non-Markovian equation 
(jHJ). Truncating the hierarchy at some term is equivalent to keeping terms 
up to a certain order in the expansion of e t ® L . In particular, if we define 
w k o(uo,t) = P f e^~ s ^ L PLe s ® L QLuo k ds, (where we have denoted explicitly 
the dependence of the quantity w k Q on the resolved initial conditions uq), 
we get 



dw k0 _ 
dt 

We can define 



Pe tL PLQLu ok + P [ e sL PLe {t ' s)QL QLQLu ok ds. 
Jo 



w k 

then 



i(u ,t) = P [ e sL PLe {t - s)QL QLQLu ok ds, 
Jo 



dw k 
dt 



rt 

- = Pe tL PLQLQLu ok + P e sL PLe {t - s)QL QLQLQLu ok ds. 

Jo 



Similarly, we can define terms w k 2, ■ ■ ■ , w kn , .... For the nth quantity w kn 
we have 



dw kn 
dt 



Pe tL PL(QL) n QLu ok + P [ e sL PLe {t ~ s)QL (QL) n+1 QLu ok ds. 

Jo 
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Of course, we have to truncate the hierarchy at some point. This amounts 
to setting the integral term appearing in the equation of evolution for w kno 
(for some index no) equal to zero. Gathering the equations up to term no 
we get 

Pe tL PLu ok + w ko (t), 
Pe tL PLQLu ok + w kl (t), 
Pe tL PLQLQLu ok + w k2 (t), 

=i = Pe tL PL(QL) n °QLu ok 

By commuting the projection P with the nonlinear functions appearing on 
the RHS of the equations, the above system of equations can be transformed 
into a closed system for Pe tL UQ k . 

4.2 Higher order models 

Higher than the second order terms in the expansion of the orthogonal dy- 
namics operator are hard to derive without some kind of book keeping rules, 
which can help reduce the risk of a mistake. There are three determining 
factors in the appearance of the expressions for a term of some order. 

1. The form of the nonlinearity, 

2. The specific form of the projection, 

3. The general property of any projection PQ = P{I — P) = 0. 

We proceed to analyze how the three factors mentioned above determine the 
form of the expressions appearing in higher order terms. This will lead us 
eventually to the formulation of a general scheme for computing higher order 
terms. First, the quadratic nonlinearity in the Euler equations and the fact 
that each new term in the series invloves one more differentiation (the L in 
the factor QL), result in the term of order n involving n + 3 powers of the 
Fourier modes. For the zeroth, first and second order models this gave cubic, 
quartic and quintic powers as already pointed out. Also, the rath order term 
will involve all combinations of n + 3 of the form n + 3 = m + I, where m, I 
are positive integers. For example, the second order term (n = 2) involves 



Jt Pe^u ok 



dw k 



o 



dt 
dw k \ 

dt 
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quintic expressions (n + 3 = 5) of the form 4 + 1, 1 + 4, 3+2, 2 + 3. In addition, 
the cubic expressions can be further decomposed as 2 + 1, 1 + 2, (see B^{uq)) 
and, of course, the trivial decomposition 3 = 3 + = + 3 which is nothing 
but Z°(uq). On the other hand, the quartic expression does not appear in 
a decomposed form but only as Z 1 (-u ). So, there is something more to 
the way the higher order terms than just a combination of all the possible 
arithmetic combinations of n + 3. In fact, there is much more structure 
that comes from the next two factors above, namely the special properties 
of the projection used and also the general property of every projection 
operator. However, before we go into more details let us make a comment: 
the structure of the expressions is convenient for numerical implementation 
because the different expressions are convolution sums in Fourier space, 
which can be computed through the Fast Fourier Transform (FFT) in real 
space, by FFT of appropriate arrays. The adjective appropriate is the catch. 
Even though all the expressions are convolution sums, one has to find what 
are the ranges for the wavenumber indices appearing in the sums. Factors 
2) and 3) will help us determine what is this range. 

To reveal this structure, we should examine more closely how a term 
in the series is related to the one preceding it. The nth order term has 
the form Z u (uq) = P L{QL) n QLuq^-, he. n applications of the operator 
QL and then application of the operator PL. We have PL(QL) n QLuok = 
PLQL{QL) n ~ l QLuQk- The part {QL) n ~ l QLuo) i is common with the n — 1st 
term Z n_1 (^o) = PL(QL) n_1 QL«o fc . When we act on {QL) n - l QLu ok with 
the extra factor QL = L—PL, we get QL(QL) n ~ 1 QLuQk = L(QL) n ~ 1 QLuok- 
Z n ~ l {uo). This is an important point. It tells us that the expression for 
QL{QL) n ^ 1 QLuQk contains 3 types of terms: i) Terms that could not ap- 
pear in Z n ~ l because of the special property of the projection which sets 
to zero expressions linear in uq^ for k G G; ii) Terms that could not appear 
in Z n ~ x due to the general property of any projection that PQ = and 
iii) Terms of the form h(uo) — (Ph)(uo), where (Ph)(uo) is any expression 
appearing in the term Z n_1 . So, we can assemble the expressions appearing 
in QL(QL) n ~ 1 QLuok into three groups according to i),ii) and iii). Then we 
can apply the operator PL once and we are done. Why is this grouping of 
expressions helpful? 

Let us start with the expressions iii). For those expressions, the final 
application of the operator PL needed to complete the calculation of Z n is 
trivial. For example, consider the case of the expression 




p+q=k 
pGG, q&F 



p+q=k 
pGG, q€F 
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which is of type 3 + 1 and appears in the term Z x (we have omitted the — i 
factor). This will give rise to a type iii) expression of the form 

^ k ■ LQLu 0p (u )A k u 0q - ^ A; • PLQLu 0p (u )A k u 0q = 

p+q=k p+q=k 

peG, qeF peG, qeF 

^ k ■ QLQLu 0p (u )A k u 0q . 

p+q=k 
peG, qeF 

Application of PL on this is trivial, since only the expression resulting from 
the application of L on QLQLuq p {uq) can survive the action of P. This is 
because the expression which comes from acting with L on uo q , i.e. 

^ k ■ QLQLuo p (u )A k R q (u ) 

p+q=k 

peG, qeF 

will be killed by the action of P, since PQ = 0. But then what remains from 
the action of PL is 

^ k ■ PLQLQLuo p (u )A k u 0q = ^ k ■ Zp(u )A k u 0q , 

p+q=k p+q=k 

peG, qeF peG, qeF 

a 4 + 1 term. 

For the types i) and ii) there is a little more work to do. We begin with 

^ k ■ LR p (u )A k u q 

p+q=k 

peFuG, qeG 

which is a 3 + 1 type i) term that appears in QLQLQLuQ k during the 
construction of Z 2 . This term could not have appeared in the term Z 1 due 
to the special property of our projection P which sets to zero UQ q with q € G. 
This term arises from the application of L on a 2 + 1 term. It is obvious 
that when acted upon with PL it will give rise to a 3 + 2 term only, since 
the associated 4+1 term would vanish, again by the special property of the 
projection. Thus, it will contribute the terms 

k ■ (Z°(u ) + B p (u ))A k R q (u ) = 

p+q=k 
peFUG, qeG 

Y k ■ Z^{uo)A k R q {u ) + Y k ■ B p (u )A k R q (u ) 

p+q=k p+q=k 

peFuG, qeG peFuG, qeG 
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where we have used the fact that (PLR p )(uq) = Z p (uq) + B p (uq). This is 
nothing more but a result of the possible decompositions of the number 3 as 
3 + or 2 + 1, i.e. Z p (ito) and B p (uq). So, we see that the 3 + 1 term which 
came from a 2 + 1 term and the special property of the projection, gives rise 
to a 3 + 2 term. 

Let us examine a type ii) term. Consider the term 

^ k- (Rp- R p )(u )A k R q (u ) = ^2 k ■ QLu 0p A k R Q (u ) 

p+q=k p+q=k 
peG, qeF peG, qeF 

which is a 2 + 2 type ii) term that appears in QLQLQLuQ k during the 
construction of Z 2 . This term could not have appeared in Z 1 because of 
the general property of any projection P that PQ = 0. It came from the 
application of L on a 2 + 1 term. It is obvious that when acted upon with 
PL, the 2 + 2 term will give rise to a 3 + 2 term, since the associated 2 + 3 
term would vanish, again by the property PQ = 0. Thus, it will contribute 
the term 

k ■ PLQLu 0p A k R q (u ) = ^ k ■ Z^(u )A k R q (u ) 

p+q=k p+q=k 

peG, qeF peG, qeF 

In summary, we can write down the following rules for the evaluation 
of the nth order term Z n (u ) = PL(QL) n QLuQ k in the Taylor series of the 
orthogonal dynamics operator: 

1. Write down the expression (QL) n ~ 1 QLuQ k . 

2. Apply QL and assemble the terms in the expression QL(QL) n ~ 1 QLuQ k 
in 3 groups: i) Terms that could not appear in Z n ~ x because of the 
special property of the projection which sets to zero expressions linear 
in liofc for k € G; ii) Terms that could not appear in Z n ~ l due to the 
general property of any projection that PQ = and iii) Terms of the 
form h(uo) — (Ph)(uo), where {PK){uq) is any expression appearing in 
the term Z n ~ l . 

3. Apply the operator PL to the type i) terms. An m + 1 term in the 
expression QL(QL) n ~ 1 QLu ok arose from an (m — 1) + 1 term and will 
give rise to an m + 2 term. The symmetric term 2 + m should also 
appear (the symmetric terms appear due to the rule of differentiating 
a product.) 
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4. Apply the operator PL to the type ii) terms. An m + I term in the 
expression QL(QL) n ~ 1 QLuQk arose from an m + (I — 1) term where 
the m term is of the form (Qh)(uo) for some function h(uo). It will 
give rise to an (m + l) + l term. The symmetric l + (m + l) term should 
also appear. 

5. Apply the operator PL to the type hi) terms. There are two cases: a) 
An (n — 2) + 1 term will give rise to an (n — 1) + 1 term with the n — 1st 
part being equal to Z n . The symmetric term 1 + (n — 1) should also 
appear; b) An m + 1 term with I ^ 1 will give rise to an (m + l) + l and 
an m+ (I + 1) term. The symmetric terms I + (m + 1) and (Z + 1) + m 
should also appear. 

6. Make sure that in the final expressions all possible decompositions of 
n + 3 into sums of two positive integers appear. All the expressions for 
the nth term in the series should be n + 3 powers of Fourier modes. 

7. As a last resort, forget about the rules and proceed with straightfor- 
ward differentiation. 

We return to the expressions (|19|) -([24 j) for Z 2 . Based on the rules above 
and the form of Z , one can see that the expressions p9 [l -(|20 )l arose from type 
i) terms, the expression (|21[) arose from a type ii) term and the expressions 
1)22(1 - 1)24(1 arose from type iii) terms. It is evident that the proliferation of 
expressions is rather rapid and, even with the rules, the amount of work to 
derive high order terms can quickly become significant. 

Note that the expansion of the orthogonal dynamics operator e s ® L and, 
more generally, of the expression PLe s ® L in a Taylor series is equivalent to 
expanding the response function of the orthogonal dynamics to the "field" 
created by the resolved modes. Thus, the expressions for the terms in the 
Taylor series of the PLe s ® L are what is known in the physics literature as 
sum rules [2] • 

One final comment about the different models. By construction, all the 
reduced models are incompressible. This is because the terms appearing in 
the models involve the incompressibility projection operator A^. 

5 Numerical results for the Taylor- Green problem 

In this section we present numerical results of the application of the different 
models to the 3D Euler equations using the Taylor-Green vortex as an initial 
condition. 
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5.1 Zeroth order model 



Note that due to the quadratic nonlinearity, the equation for each Fourier 
mode in the Euler equations contains interactions with Fourier modes of 
at most double the wavevector. This means that if one models the Fourier 
modes with wavevectors at most double than the resolved modes, then one 
obtains a closed system of equations for the resolved modes. This means 
that the ratio of the number of modes included in the set G over those in the 
set F should be at least 1. Of course, one can have a larger set G, depending 
on the computational power at hand. The more modes included in G, the 
better the modeling of the unresolved Fourier modes. The larger the ratio of 
the number of modes in G over those in F, the better the reduced model for 
the modes in F should become. We present results only for the case when 
the ratio is 1. 

The different terms appearing in the RHS of the equations for the re- 
duced models can be computed in real space using FFTs of appropriate 
arrays. Since for a reduced model of size N in each spatial direction we in- 
clude N additional unresolved modes in each direction, the arrays involved 
in the FFTs should be of size 2N . The fact that the model terms can 
be computed using the FFT makes their numerical implementation com- 
putationally efficient. Moreover, for the zeroth order model and for some 
expressions in the higher order terms of the form 

^ ] k ■ HpA^Cq, 

p+q=k 
p£G, q&F 

the FFT calculations involved are dealiased by construction and thus no 
extra (e.g. 3/2 rule dealiasing is needed. This is straightforward to see. 
For a calculation involving N modes in each direction, i.e. N/2 positive 
and N/2 negative, we perform FFTs of size 2N, i.e. N positive and N 
negative modes. But we are interested only on the RHS for the first N/2 
modes. This means (see (Hj) that to avoid aliasing (in the zeroth order 
term) we need for the total number of modes M used to satisfy the following 
inequality: —N - N/2 > N/2 - 1 - M which yields M > 2N. But 2N is 
exactly how many modes we use in the FFTs, and thus the zeroth model 
term calculation through FFTs is dealiased by construction. 

We have shown in Section |1] that the memory term can be decomposed 
into a sum of ordinary integrals. This can make the calculation of the 
memory very efficient. Unfortunately, when implemented, the models need 
to have the range of integration for the memory term reduced from [0, t] 



18 



to [to,t], otherwise the calculation becomes unstable. The value of to is 
dependent on the model and the initial condition used. However, there 
is no tuning needed. The results become better, the longer the range of 
integration, until the value of the range that leads to instability. Thus, trial 
and error is needed not to fit the results to some prescribed curve, but just 
to find when does the calculation becomes unstable. A Taylor series around 
the current instant cannot be expected to be accurate for long times in the 
past and this is the reason for the need to truncate the memory term's range 
of integration. Even though we have to truncate the range of integration, we 
can still salvage some of the efficiency gained by the transformation of the 
convolution integral into a sum of ordinary integrals. Let us assume that we 
want to calculate the integral 



where f(s) is any of the integrands appearing in the different order models' 
memory terms. Decompose the integral as 



The first integral is already computed at the last step. The second integral 
is the contribution of the current step. And the third integral is the con- 
tribution that needs to be subtracted due to the fact that the memory has 
a truncated range of integration. Note that if the integral is not truncated, 
i.e. to = tj, + At then the above decomposition is equal to 



This shows that if there is no truncation, then the computation of the mem- 
ory term is very efficient, since one needs only to add the contribution to 
the memory from the current step. We see that the truncation of the in- 
tegral forces us to keep track of the values of the integrand at the instants 
tk + At — to and tfc — to- But this means, that as the calculation progresses, 
we need to keep an array of length [to/ At], where [] stands for integer part. 
This array needs to be updated at the end of every step so that it always 
keeps the values of the integrand for the last [to/ At] steps. At every step 
we will use only two values from the array, but we need to keep the whole 
history of length [to/ At] because, the range of integration extends only for 
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to to the past. Even though the integral does not need to be computed from 
scratch at every step, as for the case of a convolution integral, the process of 
updating the array of past integrand values can be expensive for 3D calcu- 
lations. Also, the size of the array is considerable. For example, if we keep 
N modes in each direction, the size of the array is 3iV 3 x [t / At] (the factor 
of 3 comes from the 3 components of the velocity). Still, the computational 
time is half the one required if we keep the convolution integral form for the 
memory term. 

In order to study the asymptotic decay rate of the energy in the resolved 
modes, one has to evolve the system for long times. We evolved each case 
up to time t = 100, so that we have enough points to perform an accurate 
estimate of the decay rate exponent. The equations of motion for the Fourier 
modes were solved by the modified Euler method ^Zj . The integrals for the 
memory term were computed with the trapezoidal rule. The stepsize was set 
to At = 10~ 3 . We also performed an experiment with the Runge-Kutta 4th 
order method (again At = 10 -3 ) and Simpson's rule for the evaluation of 
the integrals. The results were practically the same as with the lower order 
method (the difference between the results was not larger than 10 -9 ). So, we 
decided to perform the rest of the experiments with the lower order method. 
As we mentioned before, the need to keep the values of the integrands in the 
past can be expensive. For our choice of stepsize and choice of integration 
interval for the memory term, we could only afford to use 8 3 resolved modes 
when implementing the model on a single processor workstation. 

It is important that the solutions respect the incompressibility condi- 
tion. We know that the models are incompressible by construction. The 
calculation of the integral terms involves the summation of many instants 
of the solution. Even though for short times, the solution at each instant is 
incompressible to double precision, eventually, the inaccuracy in the summa- 
tion involved in the integral terms can cause the solution to start deviating 
from incompressibility. However, for our experiments, even after 10 5 steps 
the divergence never became larger than 10 -14 . One can try to enforce the 
incompressibility condition by performing a projection on divergence-free 
fields at the end of each step. We tried that and the results did not change 
from the case without such a projection. This means that the violation of 
the divergence-free condition observed in our experiments is harmless. 

We present in Figure ^a) the evolution of the energy in the resolved 
modes E = \ Yl \uk\ 2 for the zeroth order model with 8 3 resolved modes. 

k£F 

The slope (in log-log coordinates) is a = —1.5066 ± 0.0005, which means 
that the energy decays as t a . In Figure ^b) we present results for the rate 
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(a) 



(b) 



Figure 1: (a) Energy evolution for the zeroth order model with N = 8 3 
modes, (b) Evolution of the energy decay rate. 

of energy decay dE/dt for the zeroth order model with 8 3 resolved modes. 
We set the parameter to to the value 2. This is the longest history we are 
allowed to keep to avoid an instability. 

The energy evolution shows an interesting trend. The energy decay 
seems to be organized in "waves" of activity, alternating periods of fast and 
slow decay. This behavior is absent in the t-model for the current resolution. 
Recall that the t-model is also a zeroth order model, albeit for the expansion 
of the whole memory integrand and not of the orthogonal dynamics operator 
only as the models here. In that paper, we found that one has to increase 
the number of resolved modes before the trend appears (more about the 
comparison between the two methods below). 

This organization of the energy decay is reminiscent of the phenomenon 
of intermittency, i.e. bursts of activity followed by intervals of relative in- 
action on the part of the flow. Of course, the phenomenon of intermittency 
is not only of temporal nature, but has a spatial manifestation too. This 
is exhibited as concentration of the highest vorticity in small regions of the 
flow. The trend we observe in the decay of the energy seems to assign a 
specific purpose to the vorticity. Starting from a smooth initial condition, 
we have a steepening of the gradients in the field. This means that smaller 
scales are excited until the vorticity producing mechanism runs out of steam. 
Then we enter a period of relative inaction, until there is a restart of the 
mechanism of steepening. Energy is transferred again to the smaller scales 
and so forth. This scenario continues until there is no energy left in the 
large scales. After that, the flow just disintegrates and eventually comes to 
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(a) 



(b) 



Figure 2: (a) Energy evolution for the first order model with N = 8 3 modes, 
(b) Evolution of the energy decay rate. 

a halt. The purpose of vorticity mentioned above is to regulate the transfer 
of energy to the small scales. This is reminiscent of the picture suggested 
by Moffatt, Kida and Okhitani ^5] of the vortex structures acting as the 
"sinews of turbulence" . 

5.2 First order model 

We continue our presentation of numerical results with the first order model. 
Figures [2ta) and (b) show the evolution of the energy and of the energy de- 
cay rate for a reduced system with 8 3 modes. The details of the numerical 
implementation are the same as in the case of the zeroth order model, includ- 
ing the truncation to which is again set to to = 2. As discussed in Section I57T1 
not all expressions for the first order model are dealiased by construction. 
The necessary dealiasing is done through the 3/2 rule. 

Inspection of the energy decay rate for the first and zeroth order model 
shows that the wavy structure of the energy decay has a shorter period for 
the first order model. This reflects in the energy decay rate as an increase in 
the number of spikes. Also, the slope of the energy decay (in log-log coordi- 
nates) has dropped to (3 = — 1.1379±0.0004, so the energy decays as t' 3 . The 
theoretical estimate for isotropic decaying turbulence in the limit of infinite 
Reynolds number [2D is t _1 for the case of complete self-preservation. The 
notion of complete self-preservation means that the solution is self-similar 
across all scales, from zero to infinity. For solutions that are only partially 
self-preserved, the estimate is t 7 with 7 < —1. For periodic solutions in a box 
of finite size we cannot satisfy the complete self-preservation requirement, 
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(a) 



(b) 



Figure 3: (a) Comparison of the energy evolution for the zeroth and first 
order models and the t-model with two different resolutions, (b) Rate of 
energy decay for t-model for reduced systems with 8 3 and 32 3 modes re- 
spectively. 

so that our estimate of a faster than i _1 decay becomes more plausible. In 
addition, the Taylor-Green vortex problem is not isotropic, and we do not 
know yet how lack of isotropy can affect the energy decay. During the sim- 
ulations the magnitude of the zeroth term is one to two orders larger than 
the first order term. 

In Figure Efa) we show how the energy decay of the zeroth and first 
order models with 8 3 modes compare to the energy decay for the t-model 
with 8 3 and 32 3 modes. We see that as the resolution in the t-model is 
increased it resembles more the results of the models in this paper. The 
slopes of energy decay for the t-model are —2.10 and —1.81 for the cases 
with 8 3 and 32 3 modes respectively. In Figure E[b) we show the significant 
difference in the rate of energy decay for the t-model as we increase the 
number of resolved modes. The energy decay rate for 32 3 modes exhibits 
organization of the energy decay in spikes, followed by periods of milder 
decay. It is important to note that the frequency of spikes in the 32 3 case 
for the t-model is much higher than any other case presented here. We 
expect that the same increase in the frequency of spikes will appear also 
for the zeroth and first model presented here when we perform experiments 
with larger models. However, the point of this paragraph is that the use of 
more sophisticated models can reveal, at a lower resolution, effects that need 
a higher resolution if a lower order model is used. Also, it shows that the 
t-model albeit seemingly crude is on the right track if one can afford large 



23 



(a) (b) 

Figure 4: (a) Energy evolution for the second order model with N = 8 3 
modes, (b) Evolution of the energy decay rate. 

calculations. Since the t-model is easier to implement than the models here, 
it appears as an interesting alternative to the more sophisticated models of 
this paper. 

5.3 Second order model 

We conclude our presentation of the results for the reduced models with 
the second order model. Figures EI a) and (b) show the evolution of the 
energy and of the energy decay rate for a reduced system with 8 3 modes. 
As discussed in Section I5.H not all expressions for the second order model 
are dealiased by construction. The necessary dealiasing is done through the 
3/2 rule. 

The results shown are for to = 1. It is obvious that the reduced system is 
unstable. The energy cannot grow to a value larger than the initial one. No 
matter how small a value we tried for to (from to = 2 down to to = 0.01), the 
instability is present. The instability should be related to a breakdown of 
the Taylor expansion. The energy decay rate is governed by the zeroth, first 
and second order terms. The magnitudes of the zeroth and second order 
terms are the main contributors, while the first order term (just as in the 
case of the first order model) is one order of magnitude smaller. 

6 Conclusions 

We have presented a collection of reduced models for the 3D Euler equations 
based on the Taylor expansion of the orthogonal dynamics operator that 
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appears in the Mori-Zwanzig formalism. The models up to second order 
were implemented and used to calculate the rate of energy decay for the 
Taylor-Green vortex problem. The results appear to be in good agreement 
with the theoretical estimates. The energy decay appears to be organized 
in "waves" of activity, i.e. alternating periods of fast and slow decay. 

We presented a set of rules that can facilitate the recursive calculation of 
higher order models. The rules are based on the observation that the form 
of the terms appearing in the reduced models is determined by the form of 
the nonlinearity, the special properties of the projection operator used and 
the general properties of any projection operator. 

The zero and first order models require a truncation of the range of 
integration for the memory term in order to produce solutions with decaying 
energy. However, there is no tuning needed. The results become better, 
the longer the range of integration, until the value of the range that leads 
to instability. Thus, trial and error is needed not to fit the results to some 
prescribed curve, but just to find when does the calculation become unstable. 

On the other hand, the second order model does not produce solutions 
with decaying energy no matter how small is the range of integration for the 
memory term. The numerical results suggest that second and higher order 
terms in the Taylor expansion may be more profitably used in the construc- 
tion of more elaborate approximations, e.g. Pade approximants. Of course, 
before one starts thinking in terms of Pade approximants, higher than second 
order terms in the Taylor expansion should be calculated and implemented. 
These higher order terms can potentially eliminate the instability that was 
found for the second order model. In addition, the numerical results suggest 
that a model including only a few low wavenumber modes should involve a 
long memory. This is in agreement with the main assumption behind the 
t-model presented in [TU] , 

The problem of constructing reduced models for the Euler equations has 
been, and still is, a great challenge for scientific computing. The models 
proposed here should be considered a first step in deriving models directly 
from the equations without ad hoc approximations. They are based on 
numerical and physical observations about the behavior of the solution. The 
terms appearing in the reduced models can be efficiently implemented by 
the use of the FFT on appropriate arrays. This makes the incorporation 
of the models in existing pseudospectral algorithms rather straightforward. 
We plan to apply the models in a parallel setting which will allow a better 
assessment of the properties of the flow field that is predicted by the models. 
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